Exploring the Activity Profile of TbrPDEB1 and hPDE4 Inhibitors Using Free Energy Perturbation

Human African trypanosomiasis (HAT) is a neglected tropical disease caused by the parasite Trypanosoma brucei (T.b.). A validated target for the treatment of HAT is the parasitic T.b. cyclic nucleotide phosphodiesterase B1 (TbrPDEB1). Although nanomolar TbrPDEB1 inhibitors have been obtained, their activity against the off-target human PDE4 (hPDE4) is likely to lead to undesirable clinical side effects, such as nausea, emesis, and immune suppression. Thus, new and more selective TbrPDEB1 inhibitors are still needed. This retrospective study evaluated the free energy perturbation (FEP+) method to predict the affinity profiles of TbrPDEB1 inhibitors against hPDE4. We demonstrate that FEP+ can be used to accurately predict the activity profiles of these homologous proteins. Moreover, we show how FEP+ can overcome challenges like protein flexibility and high sequence conservation. This also implies that the method can be applied prospectively for the lead optimization campaigns to design new and more selective TbrPDEB1 inhibitors.


Materials and methods
This work aims to validate the use of FEP+ to assess the activity profiles of multiple series of TbrPDEB1 inhibitors against the off-target hPDE4. The flexibility of the specific P-pocket in TbrPDEB1 constituted the challenge for FEP+ calculations. We considered it essential to analyze the degree of induced fit of the different ligand classes upon binding to address this issue.

Ligand Preparation
Data sets were collected from an in-house database (SI Table 1, 2, and 3). All PDE inhibitors considered in this study were built using the Maestro 19 interface 1 ; tautomer enumeration and protonation state assignment at experimental pH of 7.0 +/-2 was performed using LigPrep 2 .

Ligand docking
The ligands reported in SI Tables 1, 2, and 3 were docked into the TbrPDEB1 and hPDE4 structures using Glide 8 9 core-constrained docking with default parameters, using the X-ray ligand as reference (with the standard precision scoring function).
We then refined the alignment into the binding site of both PDE models for the subsequent FEP+ calculation using the Flexible Ligand Alignment tool 10 with the maximum common substructure algorithm.

FEP+ calculation
Free energy calculations were performed using the FEP+ method 11 with Schrödinger suite in the versions 2019-4 to 2020-1. Maps were generated with default settings (optimal topology).
The OPLS3e 12, 13 force field with custom parameters was used to predict the relative energy of binding of the ligands. A detailed description of the implementation of the methodology can be found elsewhere [14][15][16] . The missing torsions parameters of the ligands were computed using the Force Field Builder 17 tool.
Maps were generated with default settings (optimal topology). FEP+ jobs were run for 10 ns sampling time per λ-window for a total of 12 λ-windows. The output was analyzed with the FEP+ panel in Maestro. Replica exchange with solute Tempering (REST) was used as sampling method as reported in prior publications 16,18 (with the highest effective temperature of 1000 K for a typical perturbation with about 20 heavy atoms in the hot region) 16 The Bennett acceptance ratio method was used to calculate the free energy difference between neighboring λ windows. The simulations were performed at NPT ensemble (300 K, 1bar), using SPC water model.
Calculations were run on the GPU compute nodes with Dual Intel Xeon CPU E5-2620 v4 @2.10GHz (16 core total), and NVIDIA GeForce GTX1080 GPUs. For this study, we specifically ran the FEP+ simulations across 4 GPUs. The FEP+ calculation with the highest number of ligands in the map from Table 3 took almost 27 h to complete. It must be noted that the duration of FEP+ calculations is also dependent by the hardware used, the size of the system, the number of ligands, the number of lambda windows (default for the FEP+ implementation are as follows: 12 for neutral perturbations, 16 for core-hopping perturbations, and 24 for charge-hopping perturbations), and the length of the simulation. S5 Figure S1: Sequence alignment of human PDE4D, PDE4B, and parasite TbrPDEB1, colored by identity.